Phases of hot nuclear matter at subnuclear densities 
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Structure of hot dense matter at subnuclear densities is investigated by quantum molecular dy- 
namics (QMD) simulations. We analyze nucleon distributions and nuclear shapes using two-point 
correlation functions and Minkowski functionals to determine the phase-separation line and to clas- 
sify the phase of nuclear matter in terms of the nuclear structure. Obtained phase diagrams show 
that the density of the phase boundaries between the different nuclear structures decreases with 
increasing temperature due to the thermal expansion of nuclear matter region. The critical tem- 
perature for the phase separation is > 6 MeV for the proton fraction x = 0.5 and > 5 MeV for 
x = 0.3. Our result suggests the existence of "spongelike" phases with negative Euler characteristic 
in addition to the simple "pasta" phases in supernova cores until T < 3 MeV. 

PACS numbers: 21.65.+f,26.50.+x,97.60.Bw,61.20.Ja 
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I. INTRODUCTION 

In the process of the collapse-driven supernova [ljj, 
matter in the core experiences adiabatic compression: the 
central density increases from ~ 10 9 g cm" 3 at the be- 
ginning of the collapse to around the normal nuclear den- 
sity po = 0.165 fm -3 just before bounce; the temperature 
reaches ~ O(l) MeV at this point. 

At subnuclear densities, nuclear matter exhibits the 
coexistence of a liquid phase with a gas phase due to the 
internucleon interaction which has an attractive part. In 
the density region where nuclei are about to melt into 
uniform matter, it is expected that, at sufficiently low 
temperatures relevant to neutron star interiors [2j, the 
energetically favorable configuration of the mixed phase 
possesses interesting spatial structures such as rodlike 
and slablike nuclei and rodlike and spherical bubbles, 
etc., which are referred to as nuclear "pasta" HJ. 

This prediction is confirmed by several approaches as- 
suming nuclear shapes such as the liquid drop models 
and the Thomas-Fermi calculations Q, and is also 
confirmed without assuming nuclear shapes in the frame- 
work of the Thomas- Fermi approximation [9fl and of the 
quantum molecular dynamics (QMD) [l(J [H|- While nu- 
clear pasta at zero temperature is studied by several au- 
thors, pasta phases at finite temperatures relevant to su- 
pernova inner cores have not been studied yet except for 
a work by Lassaut et al. using the Thomas-Fermi approx- 
imation [12j and brief estimates of thermal fluctuations 
of the long- wave-length mode It is noted that, at 

temperatures of several MeV, effects of thermal fluctua- 
tions on nucleon distribution would be quite significant 
at subnuclear densities. However, the mean-field approx- 
imation such as the Thomas-Fermi [l3| and Hartree-Fock 
[lij approximation is not suitable to incorporate thermal 
fluctuations. 

Finite temperature effects lead to evaporation of nu- 



cleons from nuclear liquid region and smoothed nucleon 
density profiles. At lower temperatures where each 
nuclei fluctuates a little around an average species, a 
compressible liquid-drop model 15] with incorporating 
the temperature dependence of its bulk, surface, and 
Coulomb+lattice components provides a useful way to 
investigate the pasta phases at finite temperatures. As 
for the bulk component [ill ], binding energy, saturation 
density, and incompressibility [bn] , which are parameters 
characterizing saturation properties, decrease with in- 
creasing temperature while the temperature dependence 
of the symmetry energy is not significant f 18] . For sur- 
face component, thermal broadening of the nucleon den- 
sity profile reduces nuclear surface tension [l9l ]. Lattice 
energy is also modified by translational motion of nuclei 
[l5l |. However, at higher temperatures where the fluctua- 
tion of nuclear shape is significant, the above liquid-drop 
picture no longer holds; we have to call on some ab-initio 
method which does not assume nuclear shape. More in- 
terestingly, it might be possible that, at these temper- 
atures, the shape of the nuclear surface fluctuates and 
nuclei of various sizes and shapes coexist like colloid due 
to the entropy effect. 

These finite temperature effects can be well described 
by the methods of molecular dynamics (MD) for nucleon 
many-body systems (see, e.g., Ref. [20( for review). QMD 
[2ll ] , which is one of them, enables us to treat much larger 
systems than the other methods of MD do. Furthermore, 
at temperatures of several MeV, shell effects, which can- 
not be incorporated by QMD, are less important because 
they washed away by thermal fluctuations above ~ 3 
MeV. Thus QMD is an efficient and trustable method for 
studying nuclear matter at finite temperatures [22l . |23| . 

Pasta phases in supernova matter (SNM) are expected 
to affect the neutrino transport and hydrodynamics in 
supernova cores. Let us first note that the neutrino 
wavelengths, typically of order 20 fm, are comparable 
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to or even greater than the internuclear spacing, lead- 
ing to diffractive effects on the neutrino elastic scatter- 
ing off such a periodic spatial structure of nuclear matter 
Q. These effects, induced by the internuclear Coulombic 
correlations, would reduce the scattering rates and hence 
the lepton fraction Yl. For the bcc lattice of spheri- 
cal nuclei, such a reduction was examined by Horowitz 
p4| by calculating the associated static structure factor. 
It is also noteworthy that nonspherical nuclei and bub- 
bles are elongated in specific direction. In such direction, 
the neutrino scattering processes are no longer coherent, 
in contrast to the case of roughly spherical nuclei whose 
finitcness in any direction yields constructive interference 
in the scattering, which leads to the neutrino-trapping 
[25L |26| . The final point to be mentioned is that the 
changes in the nuclear shape are accompanied by discon- 
tinuities in the adiabatic index, denoting how hard the 
equation of state of the material is. These discontinuities 
may influence the core hydrodynamics during the initial 
phase of the collapse [l2| ■ 

In the present paper, we study the structure of hot 
dense matter at subnuclear densities within the frame- 
work of QMD. Simulations of nuclear matter with proton 
fraction x — 0.3 in addition to symmetric nuclear matter 
(x = 0.5) are performed because the typical value of the 
proton fraction for supernova matter is around 0.3 due to 
the trapping of the electron neutrinos (2|| [26| . We draw 
phase diagrams for x = 0.5 and 0.3 in the density versus 
temperature plane, which show the qualitative feature of 
the finite temperature effects on the structure of nuclear 
matter. The results of the present study would be help- 
ful to understand the real situation of the interior of the 
collapsing cores. 

The plan of this paper is as follows. In Section [Til 
we briefly explain the QMD model used in the present 
study and then discuss an effective temperature. In this 



section, a thermostatting method used in the simulations 
is also explained. In Section IIII1 we show snapshots of 
some typical nucleon distributions to discuss qualitative 
features of finite temperature effect. After that, we an- 
alyze the structure of matter using two-point correlation 
functions and Minkowski functionals, and finally, resul- 
tant phase diagrams are shown. Summary and conclusion 
are presented in Section IIV1 



II. FORMULATION 

A. Model Hamiltonian 

Simulating nuclear matter at subnuclear densities 
within the framework of QMD, we use a QMD model 
Hamiltonian developed by Maruyama et al. [22j . which 
is constructed so as to reproduce bulk properties of nu- 
clear matter and properties of finite nuclei. This model 
Hamiltonian, which describes interactions between nucle- 
ons, consists of the following six terms 

T~L = T + Vp au li + Vskyrmc + ^sym + VmD + Vcoulomb , (1) 

where T is the kinetic energy, Vp au i; is the Pauli potential 
introduced to reproduce the Pauli principle effectively, 
Vskyrme is the Skyrme potential which consists of an at- 
tractive two-body term and a repulsive three-body term, 
V sym is the symmetry potential, Vmd is the momentum- 
dependent potential introduced as two Fock terms of 
the Yukawa interaction and Vcoulomb is the Coulomb en- 
ergy including the constant contribution V p - e due to the 
Coulomb interaction between protons and electrons [28[ . 
The expressions of these terms are given as 
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where p^- means the overlap between the single-nucleon 
densities, pi(r) and Pj(r), for i-th and j-th nucleons given 
as 



Pij = / d3r P*( v ) Pj( r ) 



(8) 



(Tj is the nucleon spin and is the isospin 
(Tj = 1/2 for protons and —1/2 for neutrons) and 

Cp, go, 2>0j a, r, C s , c£\ Ci x \ /xi, /x 2 an d L are 
model parameters determined to reproduce the proper- 
ties of the ground states of the finite nuclei, especially 
heavier ones, and the saturation properties of nuclear 
matter (22|. A parameter set used in this work is shown 
in Table |U The single-nucleon densities Pi(r) and p~i(r) 
are given as 
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The quantity Pi(r) is introduced in the three body term 
of Skyrme interaction Eq. (j4]) to incorporate the effect 
of the repulsive density-dependent term by the modified 
width L. 



TABLE I: Effective interaction parameter set 
(/Y=280 MeV; medium EOS model in Ref. [H]) 



a (MeV) 


-92.86 


P (MeV) 


169.28 


T 


1.33333 


Cs (MeV) 


25.0 


Cfi 5 (MeV) 


-258.54 


CS 5 (MeV) 


375.6 


Mi (MeV) 


2.35 


/x 2 (MeV) 


0.4 


L (fm 2 ) 


2.1 



B. Effective temperature 

The effective Hamiltonian (TTJ used in this work con- 
tains momentum dependent interactions, i.e. the Pauli 
potential Vp aun and the momentum dependent potential 



Vmd ■ Thus the usual expression for the instantaneous 
kinetic temperature Tjdn given as 
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loses its meaning of the temperature in thermodynamics. 
We here use an effective temperature T e g proposed by 
Chikazumi et al. [29j |. which is given as 



N ^ 2 l ' dt 
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It can be immediately seen that this expression is equiva- 
lent to the usual kinetic temperature defined by Eq. (fT"2"]) 
if the effective Hamiltonian % does not have momentum 
dependent interactions (i.e., Vp au ii and Vmd)- 

In order to confirm whether the effective temperature 
T e g is consistent with temperature in the Boltzmann 
statistics , we pe rform Metropolis Monte Carlo (MC) sim- 
ulations |30l - l32l ] with 256 nucleons (with x = 0.5, i.e., 128 
protons and 128 neutrons). We investigate six given tem- 
peratures T sct =0.1, 1, 3, 5, 7 and 10 MeV at four nu- 
cleon densities p = 0.1, 0.3, 0.5 and 0.8 po within a wide 
region of the phase diagram at subnuclear densities cov- 
ering from a phase-separating region to a uniform fluid 
region. We prepare a cubic box, which is imposed of the 
periodic boundary condition. In the simulations the sys- 
tem is equilibrated at a given temperature T set for 1000 
MC steps (i.e. 1000 x N trial moves), and then sampling 
is carried out for the following 10000 MC steps. Sampled 
values of the instantaneous effective temperature T e g is 
plotted in Fig. [TJ We can see from this figure that the in- 
stantaneous effective temperature T e ff fluctuates around 
the given value of T sct . It is noted that the long-time 
averaged values (T e g) of the effective temperature coin- 
cide with r set quite well within the fluctuations of order 
~ T sct /VN due to the finite particle number. Thus we 
can conclude that the effective temperature given by Eq. 
(|13[) is consistent with temperature in the Boltzmann 
statistics. It is also confirmed that, in microcanonical 
molecular dynamics simulations, the mean value of the 
effective temperature keeps constant after the system is 
equilibrated enough. 

The instantaneous effective temperature can be neg- 
ative as plotted in Fig. Q] for T sct = 0.1 when Pi and 
Ri take the opposite directions each other due to the 
contribution of the momentum dependent interactions. 
However, it is also confirmed that, after the system is 
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MC step 

FIG. 1: Sampled values of the instantaneous effective tem- 
perature Tea and their long-time averages (T e ff) for p — 
0.1, 0.3, 0.5 and 0.8 p . 



relaxed, the long-time average of the effective tempera- 
ture does not take negative values when we pursue the 
time evolution of the system by the QMD equations of 
motion even though some friction terms were attached to 
them like Eqs. (17) in Ref. [11]. In the remaining part 
of this paper, we measure temperature T by the effective 
temperature. 



C. Thermostat 



the extended system in this case can be written as 

= £S + " ({R ' UP ' l) + W +9 T 

where U({Ri},{Pi}) is the potential which depends on 
both positions and momenta, s is the additional dynam- 
ical variable for time scaling, p s is the momentum con- 
jugate to s, Q is the thermal inertial parameter corre- 
sponding to a coupling constant between the system and 
thermostat [in our simulations, we set Q ~ 10 8 MeV 
(fm/c) 2 ], g is a parameter to be determined as 3N by 
a condition for generating the canonical ensemble in the 
classical molecular dynamics simulations, and /? is de- 
fined as f3 = 1/fceTset- The equations of motion yield 

dt ~ dP t ~ mi dPi' [ ' 



with 




where £ means the thermodynamic friction coefficient. 

In the dynamical process described by these equations, 
'Hnosc is conserved and the value of the effective temper- 
ature fluctuates around T SGt as can be seen from Eq. (|18p . 



It is necessary to perform "isothermal" QMD simu- 
lations in order to equilibrate the system at a specified 
effective temperature. In ordinary molecular dynamics 
simulations, what is called the Nose-Hoover thermostat is 
commonly used to carry out constant-temperature sim- 
ulations |3lT[34| . The approach of Nose and Hoover is 
based on the Hamiltonian of an extended system, which 
contains additional and artificial coordinates and veloc- 
ities intended to mimic the dynamics of the system in 
contact with a thermal bath. It is shown by Nose that 
this method generates the states in the canonical ensem- 
ble average, i.e., the microcanonical ensemble average in 
the extended system reduces to the canonical ensemble 
average in the real system. For example, the momen- 
tum distribution coincides with the Maxwell-Boltzmann 
distribution exactly after the system is well equilibrated 
with this thermostat. 

Here we modify the Nose-Hoover's method so as to 
adapt to the effective temperature. The Hamiltonian of 



III. SIMULATIONS AND RESULTS 
A. Procedures for simulations 

Let us here explain the procedure for the simula- 
tions. We investigate nuclear matter with proton frac- 
tions x = 0.5 and 0.3 at subnuclear densities in suffi- 
ciently wide regions of the density versus temperature 
plane covering the whole region where phase separation 
is observed: symmetric nuclear matter is studied up to 
p = 0.7po and T = 8 MeV, and nuclear matter with 
x = 0.3, up to p = 0.6po and T = 7 MeV. Intervals 
of the density and the temperature between the inves- 
tigated points are 0.025po or 0.05po and 0.5 MeV or 1 
MeV, respectively (from the present section onward, we 
set fee = 1). 

We perform simulations for a cubic box with periodic 
boundary condition. We study the (n,p,e) system with 
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2048 nucleons (for some typical cases of the columnar 
phase and the planar phase, a system with 16384 nucle- 
ons is also used). Throughout this paper, we treat sys- 
tems which are not magnetically polarized, i.e., they con- 
tain equal numbers of protons (and neutrons) with spin 
up and spin down. The relativistic degenerate electrons 
which ensure charge neutrality are regarded as a uniform 
background because the influence of the electron screen- 
ing on the phase diagram at subnuclear densities is small 
as shown explicitly in Ref . [13] • The Coulomb interaction 
is calculated by the Ewald method taking account of the 
Gaussian charge distribution of the proton wave packets 
(see, e.g., Appendix A in Ref. [HI), which enables us to 
sum up the contributions of long-range interactions in a 
system with periodic boundary conditions efficiently. For 
the nuclear interaction, we use the effective Hamiltonian 
developed by Maruyama et al. (medium EOS model) 
22] , whose expressions are given in Section IlLAl 

We first prepare a hot, uniform gas with 2048 nucleons 
at T ~ 20 MeV as an initial condition, which is equili- 
brated for ~ 500 — 2000 fm/c in advance. We then cool 
it down slowly for O(10 3 — 10 4 ) fm/c keeping the nucleon 
density unchanged by the frictional relaxation method 
[see Eqs. (17) of Ref. [Hj |] until the temperature reaches 
~ 5 MeV. For the present QMD model, this is the typical 
temperature for the boundary of the phase-separating re- 
gion at subnuclear densities relevant to the pasta phases. 
In some cases, the thermostat of the Nose-Hoover type 
(see Section III C| is also used to cool the system down 
quickly until ~ 10 MeV, at which temperature matter is 
still completely uniform. 

After the cooling process, the system is then relaxed 
for ~ 4000 — 5000 fm/c at a given temperature T sct us- 
ing the thermostat of the Nose-Hoover type, which is 
followed by a further relaxation for ~ 5000 fm/c at the 
same T set without the thermostat (i.e., microcanonical 
molecular dynamics simulation). Thermal averages are 
measured in the microcanonical relaxation process. The 
above relaxation processes with and without the thermo- 
stat are repeated for the other values of T set by changing 
T sct by 0.5 or 1 MeV and keeping the density constant. 

Simulations of a larger system with 16384 nucleons 
have also been performed for some typical cases of the 
phases with slablike nuclei and with rodlike nuclei to ex- 
amine the importance of finite size effects. We combine 
eight replicated samples at T = with 2048 nucleons into 
a 16384-nucleon sample. We then add numerical noise to 
the positions and the momenta of nucleons, up to 0.1 fm 
in the position and 1 MeV/c in the momentum. We in- 
crease the temperature by lMeV and relax the system 
for ~ 4000 — 5000 fm/c using the Nose- Hoover thermo- 
stat and relax further for ~ 3000 — 5000 fm/c without 
the thermostat. These relaxation processes are repeated 
for T sct = 2 and 3 MeV. 

The simulations of the 2048-nucleon system are 
performed using PCs (Pentium III) equipped with 
MDGRAPE-2, and those of the 16384-nucleon system are 
done by Fujitsu VPP 5000 equipped with MDGRAPE-2. 



B. Two-point correlation functions and Minkowski 
functionals 



To analyze the spatial distribution of nucleons, we use 
the two-point correlation function. The two-point cor- 
relation function £u for the nucleon density field 
(i = N,p,n; where N stands for nucleons) is here de- 
fined as 

&( r ) = hj^vf ^^(x^x + r) (20) 
= (^(x^x + r))^ , (21) 

where (• • -) Xl n r denotes an average over the position x 
and the direction of r, and <5j(x) is the fluctuation of the 
density field p^ (x) given by 



_ pW(x)-pW 
Oi(x) = == 



with 



V 



(22) 



(23) 



To identify the nuclear surface and extract its morpho- 
logical characteristics, we use the Minkowski functionals 
(see, e.g., Ref. [35| and references therein; a concise re- 
view is provided by Ref. 36] ; a brief explanation is given 
in Section IV C of Ref. ll]), especially of the integral 
mean curvature and the Euler characteristic [53] • 

Suppose we set a threshold density pth and consider the 
regions where the density is higher than this value sur- 
rounded by the isodensity surfaces for p t h (the procedure 
for identifying the nuclear surface, which is characterized 
by isodensity surfaces for a specific value of pth, will be 
explained in Section MIDI) . The integral mean curvature 
and the Euler characteristic arc defined as surface inte- 
grals of the following local quantities: the mean curvature 
H = (ki + K2)/2 and the Gaussian curvature G = K1K2, 
i.e., J gK HdA and x = ^ Jg^ GdA, where K\ and K2 are 
the principal curvatures and dA is the area element of 
the surface of the body K. The Euler characteristic \ is 
a purely topological quantity and is expressed as 

X = (number of isolated regions) — (number of tunnels) 
+ (number of cavities). (24) 

Here we introduce their normalized quantities: the area- 
averaged mean curvature, (H) = -j J HdA, and the Eu- 
ler characteristic density, x/V, where V is the volume of 
the whole space. 

In the present work, we use 64 3 (128 3 ) grid points for 
the 2048-nucleon (16384-nucleon) system in constructing 
the nucleon density distribution pW(x). Detailed pro- 
cedures for calculating £y (( g) and x/V) are given in 
Section IV A (IV C) of Ref. [n|. 
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(a)T=1MeV (b) T=2MeV (c) T=3MeV (d) T=4MeV 




FIG. 2: (Color online) The nucleon distributions for x = 0.5, p = 0.225/9o at the temperatures of 1,2,3 and 4 MeV. 2048 
nucleons are contained in the simulation box of size Lbox = 38.07 fm. These figures show the top views along the axis of the 
rodlike nuclei. Protons are represented by the red particles, and neutrons by the green ones. 



(a) T=1 MeV (b) T=2MeV (c) T=3MeV 







FIG. 3: (Color online) The same as Fig.[2]at the temperatures of 1,2 and 3MeV for the system with 16384 nucleons. The box 
size Lbox is 76.14 fm. The upper panels show the top views along the axis of the cylindrical nuclei at T = 0, the lower ones the 
side views. 



C. Typical nucleon distributions for the phases 
with rodlike and slablike nuclei 



Let us first show some snapshots of the nucleon distri- 
bution at finite temperatures for densities corresponding 
to the phases with slablike nuclei and with rodlike nuclei 
at T — 0. These snapshots help us to understand the 
qualitative feature of finite temperature effects on the 
nuclear structure. 



Figures [5] and [3] show snapshots of the nucleon distri- 
bution for x = 0.5 at a density of 0.225po (the phase with 
cylindrical nuclei at T = 0) and Figs. 2] and [5] show those 
for x — 0.5 and p — 0.4po (the phase with planar nuclei 
at T = 0). Here we show snapshots of the 2048-nucleon 
system and those of the 16384-nucleon system for both 
cases. Figures [8] and |U for x = 0.3 are the same as Figs. 
[3] and [5] for x = 0.5, respectively; Fig. [8] is for 0.175po 
(the phase with cylindrical nuclei at T = 0), and Fig. 
[9] is for 0.34po (the phase with planar nuclei at T = 0). 
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(a)T=1MeV (b) T=2MeV (c) T=3MeV (d) T=4MeV 




FIG. 4: (Color online) The nucleon distributions for x — 0.5, p = 0.4po at the temperatures of 1,2,3 and 4 MeV. 2048 nucleons 
are contained in the simulation box of size Lbox = 31.42 fm. Protons are represented by the red particles, and neutrons by the 
green ones. 



(a) T=1 MeV (b) T=2MeV (c) T=3MeV 




FIG. 5: (Color online) The same as Fig.|4]at the temperatures of 1,2 and 3MeV for the system with 16384 nucleons. The box 
size Lbox is 62.84 fm. These figures are shown in the direction parallel to the plane of the slablike nuclei at T = 0. 



The snapshots for the 16384-nucleon systems (Figs. [3j EJ 
and El) are depicted without perspective. 

From these figures, we can see the following qualitative 
features irrespective of the proton fraction and the sys- 
tem size: at T ~ 1.5 — 2 MeV (but snapshots for T ~ 1.5 
MeV are not shown there), the number of the evapo- 
rated nucleons starts to be significant; at T > 3 MeV, 
the density profiles of the nucleons are smoothed out and 
it is difficult to identify the nuclear surface. In view of 
the fact that these general features are the same for sys- 
tems with different particle number (see Figs. [2] and [3] for 
p = 0.225po; Figs. 2] and 03 for p = OApo), we can say that 
a qualitatively correct phase diagram can be obtained by 
using 2048-nucleon system. This statement is supported 
by the behaviors of the two-point correlation function 
for these two systems. It is remarkable that £nn for the 
2048-nucleon system and that for 16384-nucleon one co- 
incide quite well, as shown in Figs. [6] and [7] However, 
we should note that the larger system can incorporate 
thermal fluctuations of longer wavelengths (38[. As can 
be seen by comparing Figs. SJa) and 0a), the slablike 
nuclei have waves in the 16384-nucleon system at T = 1 
MeV, but they do not in that with 2048 nucleons. 



By comparing the two cases of different values of the 
proton fraction in more detail, we can see that the num- 
ber of the evaporated protons is significantly smaller than 
that of neutrons at T ~ 2McV for x — 0.3 although they 
are close to each other at T ~ 2MeV for x = 0.5. Here, 
the number density of the evaporated neutrons is defined 
as the number density of neutrons outside nuclei minus 
that of dripped neutrons at T = 0. A nuclear matter 
region with higher proton fraction (but less than 0.5) is 
more energetically favorable than that with lower pro- 
ton fraction because of the symmetry energy, and thus, 
at x = 0.3, neutrons are preferentially evaporated to in- 
crease the proton fraction in the nuclei. 

We also note that, for x = 0.5, the slablike nuclei touch 
and fuse with each other at T = 2 MeV [see Figs. 0Jb) 
and0b)] while the rodlike structure persists at this tem- 
perature [see Figs. EJb) andHJb)]. The fragility of the 
phase with slablike nuclei would stem from the Landau- 
Peierls instability and its larger volume fraction for the 
nuclear matter region. 

In closing the present section, we would like to mention 
the effect of the dripped neutrons on the nuclear struc- 
ture. Here we note that the dripped neutrons, on the one 
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hand, suppress the thermal expansion of nuclei due to its 
pressure acting on the nuclear surface, and they, on the 
other hand, reduce the nuclear surface tension. As can 
be seen by comparing Figs. Ela), (b) andUJa), (b), rod- 
like nuclei tend to buckle at x = 0.3 as the temperature 
increases, but they do not at x = 0.5 while just expand 
in radius. At a fixed density and a fixed number of nu- 
clei, the nuclear radius is directly related to the volume 
fraction of the nuclear matter region, which is principally 
determined by the bulk properties. Conversely, bending 
of nuclei is controled by elastic constants of nuclei, which 
depend on nuclear surface tension and the Coulomb en- 
ergy [3!| . The reduction of nuclear surface tension due 
to the dripped neutrons leads to decrease of the elastic 
constants, which makes rodlike nuclei easy to buckle in 
the case of x = 0.3. In the case of x = 0.5, there is no 
reduction of the elastic constants due to the dripped neu- 
trons, and we thus only observe the thermal expansion in 
radius of the rodlike nuclei. Next, let us consider about 
slablike nuclei. By comparing Figs. [5jb) and^b), we 
can see that, for x — 0.5, the slablike nuclei expand in 
thickness and they touch with each other at T ~ 2 MeV; 
for x = 0.3, their expansion is smaller and they do not 
touch at this temperature. This result can be understood 
as a consequence of suppressing the thermal expansion of 
nuclei due to the dripped neutrons. 
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FIG. 7: The same as Fig. [6] for x = 0.5 and p = 0.4p o , 
where the system is in the phase with rodlike nuclei at zero 
temperature. 



D. Phase diagrams 




r (fm) 

FIG. 6: Comparison of the two-point correlation functions 
£nn for the systems with 2048 and 16384 nucleons. These are 
calculated for x = 0.5 and p = 0.225po, where the system is in 
the phase with rodlike nuclei at zero temperature. The error 
bars are the standard deviations in the long-time average. 



In constructing phase diagrams, we determine the 
phase separating region, identify the nuclear surface and 
classify the nuclear structure according to its morpholog- 
ical characteristics. For these purpose, we use the two- 
point correlation functions and the Minkowski function- 
als introduced in Section UlIBI 

The phase-separation region is determined by the 
two-point correlation function £nn of the nucleon den- 
sity fluctuation. In Figs. [TOl and [Til we plot £jVJv( r ) 
for symmetric nuclear matter at p — 0.225 and 
OApo, respectively, as typical examples. The former den- 
sity corresponds to the phase with rodlike nuclei and the 
latter one to that with slablike nuclei at zero temperature 
(see Figs. [5] and HJ). From these figures, we can see that 
the amplitude of £jv n an d the relative density dispersion 
given by £nn(0) decrease as the temperature increases. 
It is also noted that the smallest value of r at which £jvjv 
is zero does not change much compared to the change 
in the amplitude. This behavior of £,nn is similar to 
that when the density increases keeping T = 0, which we 
have studied in detail in Ref. fTlj ] . The important point 
is that, from the behavior of ^nn, we can determine the 
temperature at which the long-range correlation of the 
nucleon distributions disappears. Figure [TU] (jlip shows 
that the long-range correlation cannot be seen at T = 6 
MeV (5 MeV) and thus the phase-separating boundary 
lies between T = 5 and 6 MeV at 0.225p (T = 4 and 5 
MeV at 0.4p )- 

The two-point correlation functions £u (i = N,p, and 
n) for x = 0.3 are also plotted in Figs. [T2] and [JjJ] for 
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(a)T=1MeV (b) T=2MeV (c) T=3MeV 




FIG. 8: (Color online) The nucleon distributions for x = 0.3, p = 0.175po at the temperatures of 1,2, and 3 MeV. 16384 
nucleons are contained in the simulation box of size Lbox = 82.788 fm. Protons are represented by the red particles, and 
neutrons by the green ones. The upper panels show the top views along the axis of the cylindrical nuclei at T = 0, the lower 
ones the side views. 



(a) T=1 MeV (b) T=2MeV (c) T=3MeV 




FIG. 9: (Color online) The nucleon distributions for x = 0.3, p — 0.34po at the temperatures of 1,2, and 3 MeV. 16384 nucleons 
are contained in the simulation box of size Lbox = 66.34 fm. Protons are represented by the red particles, and neutrons by the 
green ones. These figures are shown in the direction parallel to the plane of the slablike nuclei at T = 0. 



p = 0.175 and 0.35/9 , respectively. As shown in Figs. [5] 
andUl these densities correspond to the phases with rod- 
like nuclei and with slablike nuclei (Although the snap- 
shots shown in these figures are for the 16384-nucleon sys- 
tem, the correlations functions £u shown in Figs. [12] and 
[T3]are calculated for the 2048- nucleon system.). The am- 
plitude of £jj and the relative density dispersion decrease 
with increasing temperature as in the case of a; = 0.5. 
Among fjvjv: Cppi an d £,nm the second one has the largest 
amplitude and the third one has the smallest except for 



high temperatures at which matter is almost uniform. 
We note that the smallest zero-points of the three kinds 
of £a almost coincide each other and they remain nearly 
constant at lower temperatures. This behavior shows 
that the density fluctuations of protons and neutrons are 
strongly correlated even at x — 0.3 and the wavelength 
of the density fluctuations does not change significantly 
at these temperatures. As can be seen from Fig. [T2Ta) 
[r3Ta)]. the long-range correlation disappears at T > 5 
MeV [T = 4 MeV] and thus the phase-separating bound- 
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FIG. 10: Two-point correlation function of the density fluc- 
tuation calculated for N = 2048, x = 0.5 and p = 0.225po, 
where the system is in the phase with rodlike nuclei at zero 
temperature. The error bars are the standard deviations in 
the long-time average. 




5 10 15 20 

r (fm) 

FIG. 12: Two-point correlation function of the density fluc- 
tuation calculated for N = 2048, x = 0.3 and p — 0.175po, 
where the system is in the phase with rodlike nuclei at zero 
temperature. The error bars are the standard deviations in 
the long-time average. 




r (fm) 



FIG. 11: Two-point correlation function of the density fluctu- 
ation calculated for N = 2048, x = 0.5 and p = 0.4po, where 
the system is in the phase with slablike nuclei at zero tem- 
perature. The error bars are the standard deviations in the 
long-time average. 




r (fm) 

FIG. 13: Two-point correlation function of the density fluc- 
tuation calculated for N — 2048, x = 0.3 and p = 0.35po, 
where the system is in the phase with slablike nuclei at zero 
temperature. The error bars are the standard deviations in 
the long-time average. 
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ary lies between T = 4 and 5 MeV at 0.175/9 [T = 3 and 
4 MeV at 0.35p ]- 
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FIG. 14: Euler characteristic density for various temperatures 
plotted as a function of the threshold density p t h- The data 
shown here are calculated for x — 0.5 and p = 0.225po, where 
the system is in the phase with rodlike nuclei at zero tem- 
perature. The error bars are the standard deviations in the 
long-time average. For T = 2 MeV, the plateau region lies 
between p t h — 0.45po and ~ 0.6po, where x/V is negative. 
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FIG. 16: Euler characteristic density for various temperatures 
plotted as a function of the threshold density p t h. The data 
shown here are calculated for x = 0.3 and p — 0.175po, where 
the system is in the phase with rodlike nuclei at zero tem- 
perature. The error bars are the standard deviations in the 
long-time average. 
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FIG. 17: The same as Fig. [16] for 0.35po, where the system is 
in the phase with slablike nuclei at zero temperature. 



FIG. 15: The same as Fig.[l4]for p — 0.4po, where the system 
is in the phase with slablike nuclei at zero temperature. 

To classify the phase in terms of the nuclear struc- 
ture, we first have to identify nuclear surface. The Euler 
characteristic is useful for identifying the nuclear surface. 
Figures [14] and Qj)] show the Euler characteristic density 
x/V of the nucleon density field calculated for symmetric 
matter for p = 0.225 and 0.4p , respectively, and they are 
plotted as a function of the threshold density p t h for an 
isodensity surface. As for x = 0.3, x/V is plotted in Figs. 
fl~6l and HTl for p = 0.175 and 0.35,oo, respectively. It is 
noted that the curves of x/V for lower temperatures have 



a plateau, which corresponds to the nuclear surface. The 
plateau value of x/V is in the range of — O(10 4 ) to O(10 4 ) 
(see, e.g., Figs. 12 and 16 in Ref. [Hj]), * nus we should 
use the region where the error is smaller than ~ 10~ 4 for 
judging whether x/V is zero or not (hereafter, we use the 
word "plateau" for the plateau region in which the error 
is smaller than ~ 10~ 4 ). According to these figures, the 
plateau of x/V can be observed at T = 1 and 2 MeV, 
while it cannot be seen for T > 3 MeV. Thus we can say 
that the nuclear surface cannot be identified in the tem- 
perature range of T > 3 MeV at these values of x and 
p even in the phase-separating region. This result agrees 
with the qualitative feature of the snapshots of nucleon 
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distributions shown in Figs [2] [5J [8] and [9l 

In the region where the nuclear surfaces can be iden- 
tified, we sort them into several phases by their nuclear 
shapes. We extract morphological features of the nu- 
clear surface using the Euler characteristic density x/V 
and the averaged mean curvature (H). As one can see 
from Eq. (l24l) . x/V > f° r the sphere and the spherical 
hole phases, which have isolated regions and cavities, re- 
spectively (the coexistence phase of spheres and cylinders 
also has x/V > 0), and x/V = for the other ideal pasta 
phases (see the plateau of x/V = f° r lower tempera- 
tures shown in Figs. HU [15j [TH and IT7| . i.e., the cylin- 
der, the slab and the cylindrical hole phases which con- 
sist of infinitely long rods, infinitely extending slabs and 
infinitely long cylindrical holes, respectively. The "spon- 
gelike" intermediate phases which are obtained in our 
previous works about zero-temperature matter fill fill ] 
have negative values of x/V (see Fig. 12 of Ref. |llj|). 

As for the averaged mean curvature, at zero tempera- 
ture, this quantity decreases almost monotonically from 
positive to negative with increasing density until matter 
becomes uniform (see Fig. 12 of Ref. [llj). Obviously, 
the zero-value of (H) corresponds to the slab phase. 

In our previous work [llj ]. the sequence of nuclear 
shapes at zero temperature with increasing density is de- 
scribed by the quantities (H) and x/V as follows: sphere 
((H) > 0, x/V > 0) -> cylinder ((H) > 0, X /V = 0) 
-> ((H) > 0, x/V < 0) -> slab ((H) = 0, x/V = 0) -> 
((H) < 0, x/V < 0) -> cylindrical hole ((H) < 0, x/V = 
0) -> spherical hole ((H) < 0, x/V > 0) — !• uniform. 

Here we classify the nuclear structure at various tem- 
peratures and densities according to the combinations of 
signs of (H) and x/V as we have done in the case of 
T = [1J|. The phase diagrams obtained for x = 0.5 
and 0.3 are plotted in Figs. [TS] and Qj|l respectively, in 
the p - T plane 40] . Qualitative features of these phase 
diagrams are the same, but the phase separating-region 
surrounded by a dashed line is smaller for x = 0.3 than 
that for x — 0.5. As shown by dotted lines in these fig- 
ures, nuclear surface can be identified typically at T < 3 
MeV in the density range of interest. In this region where 
we can observe the nuclear surface, the nuclear struc- 
ture generally changes in the same sequence mentioned 
above even at finite temperatures (the only exception is 
that the density region of the phase with slablike nuclei 
is bounded at T > 2 MeV for x = 0.3). The regions 
above the dotted line and the dashed line correspond to 
some non- uniform phase, which is however difficult to be 
classified into specific phases because the nuclear surface 
cannot be well identified. 

The most striking finite temperature effect on the 
phase diagrams is that the density of each phase bound- 
ary between the different structures decreases as T in- 
creases. This is due to the thermal expansion of the 
nuclear matter region. At higher temperatures, the vol- 
ume fraction of the nuclear matter region is larger even 
at smaller densities, which leads to decrease of the each 
phase boundary (at 0.225po for % — 0.5 and at 0.175/9q 



for x = 0.3, the volume fraction of the nuclear matter 
region increases by ~ 5% at T = 2 MeV from the value 
at T = 1 MeV; at 0Ap for x = 0.5 and at 0.35p for 
x = 0.3, this increase is ~ 10 — 15%). This feature is less 
noticeable for x = 0.3 than x = 0.5, probably because the 
thermal expansion of nuclei is prevented by the pressure 
of the dripped neutrons outside of nuclei. 

Above T ~ 3 MeV, single-particle excitations are dom- 
inant rather than collective nuclear deformation, and si- 
multaneously, the binding potential of "nuclei" becomes 
smeared out. Further increase of temperature leads to 
dissolution of nuclei into the uniform fluid phase. The 
critical point for the phase separation is located at p ~ 
0.25p and T > 6 MeV for x = 0.5 and at p ~ 0.2/9 and 
T > 5 MeV for x = 0.3. The above value of the critical 
temperature for x = 0.5 is not so different from the re- 
sult obtained by Chikazumi et al., T c — 8 MeV, using a 
similar nuclear potential [2j|. 

As far as we have investigated, the colloidal state in 
which nuclei of various sizes and shapes coexist cannot 
be observed except for the coexistence phase of spheres 
and cylinders. To resolve the problem whether the col- 
loidal phases are realized due to entropy effect at finite 
temperatures, it would be necessary to perform simula- 
tions with the 16384-nucleon system for various densities 
other than the typical cases shown in the Section llll CI or 
for the temperature region between 2 and 3 MeV care- 
fully, in which the nuclear surface still remains but its 
fluctuations are large. 

IV. SUMMARY AND CONCLUSION 

We have investigated the structure of hot nuclear mat- 
ter at subnuclear densities by QMD simulations. The val- 
ues of the proton fraction studied in the present work are 
0.5 and 0.3; the latter is a typical value in supernova in- 
ner cores. The phase diagram obtained for x — 0.3 (Fig. 
IT51) tells that the pasta phases and the spongelike inter- 
mediate phases survive until T < 3 MeV except for the 
higher density region of p > 0.45po. This result strongly 
suggests that these phases with exotic nuclear structures 
can exist in the inner cores in the pre-bounce stage of the 
collapse-driven supernovae. Thus it would be meaning- 
ful to investigate a modification to the neutrino opacity 
due to the infinitely extending "pasta" and spongelike 
structures in the future f4lj . 

The present research has provided the general picture 
of the phase diagram at subnuclear densities in the tem- 
perature vs density plane. The qualitative features of 
the nucleon distribution at finite temperatures, except 
for higher densities just below which matter becomes uni- 
form even at T ~ MeV, may be summarized as follows. 

• At T ~ 1-1.5 MeV, the number of evaporated nu- 
cleons starts to be significant. However, the macro- 
scopic structure of matter does not basically change 
from that of zero-temperature matter except for a 
small expansion of the nuclear matter regions. 
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FIG. 18: Phase diagram of matter at x = 0.5 plotted in the p - T plane. The dashed and the dotted lines on the diagram 
show the phase separation line and the limit below which the nuclear surface can be identified, respectively; thus the regions 
above the dotted line and the dashed line show the non-uniform regions, which are difficult to be classified into specific phases 
in terms of the nuclear structure. The dash-dotted lines are the phase boundaries between the different nuclear shapes. The 
symbols SP, C, S, CH, SH, U stand for nuclear shapes, i.e. sphere, cylinder, slab, cylindrical hole, spherical hole and uniform, 
respectively. The parentheses (A,B) show intermediate phases between A and B-phases. The regions (a)-(g) correspond to the 
nuclear shapes characterized by (H) and x/V as follows: (a) (H) > 0, x/V > 0; (b) (H) > 0, x/V = 0; (c) (H) > 0, x/V < 0; 
(d) (H) = 0, x/V = 0; (c) (H) < 0, x/V < 0; (f) (H) < 0, x/V = 0; (g) (H) < 0, x/V > 0. Simulations have been carried 
out at points denoted by circles. 
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• At T ~ 1.5-2.5 MeV, fluctuations of the nuclear 
shape start to be significant and, in some cases, 
the topology of the nuclear structure changes from 
that in the zero-temperature matter for the same 
density. 

• At T ~ 2.5 - 3 MeV, single particle excitations are 
dominant rather than collective excitations, repre- 
senting nuclear deformations. The density profile of 
nucleons is smoothed out, and nuclear surface can- 
not be identified by a threshold density. At higher 
temperatures, even the density contrast becomes 
drastically smaller, which leads to the transition to 
uniform matter. The order of the transition from 
inhomogeneous matter to uniform matter has yet 
to be clarified. 

• The critical temperature T c for the phase separa- 
tion is > 6 MeV for x = 0.5 and > 5 MeV for 
x = 0.3. 

Through the present work and the previous ones 
[Tol Hil ] , we have depicted the structure of nuclear matter 
in neutron star inner crusts and supernova inner cores in 
its equilibrium state. Our studies will be helpful in un- 
derstanding the realistic situations of the interior of the 
dense stars. For further work, studies on the dynamical 
aspects of the nuclear pasta become important. For ex- 



ample, in connection with the possible instability of the 
long-range order of the planar phase, which cannot be 
fully incorporated in our simulations, elastic properties 
[39| of the pasta and spongelike phases should be studied 
further. Structural transitions between the pasta phases 
induced by compression and decompression is one of the 
most interesting problems (42j which have not yet been 
studied at all. We believe that our work here, which was 
performed within a dynamical framework, will open up 
these new aspects of the study on nuclear pasta. 
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